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Abstract 

We consider a general class of high order weak approximation schemes 
for stochastic difTerential equations driven by Levy processes with infi- 
nite activity. These schemes combine a compound Poisson approximation 
for the jump part of the Levy process with a high order scheme for the 
Brownian driven component, applied between the jump times. The over- 
all approximation is analyzed using a stochastic splitting argument. The 
resulting error bound involves separate contributions of the compound 
Poisson approximation and of the discretization scheme for the Brownian 
part, and allows, on one hand, to balance the two contributions in order 
to minimize the computational time, and on the other hand, to study 
the optimal design of the approximating compound Poisson process. For 
driving processes whose Levy measure explodes near zero in a regularly 
varying way, this procedure allows to construct discretization schemes 
with arbitrary order of convergence. 
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1 Introduction 

Let Xt be the unique solution of the SDE 

Xt = x+ [ b{X,)ds+ [ cj{Xs)dBs+ I h{X,_)dZs, (1) 
Jo Jo Jo 

where b, a and h are functions with bounded derivatives, iJ is a (multi- 
dimensional) Brownian motion and Z a one-dimensional infinite activity pure 
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jump Levy process with Levy measure v. In this paper we are interested in the 
weak approximation of Xj, using random partitions of the time interval. 

The traditional approach, analysed, e.g., in Jacod et al |Tn| and Protter- 
Talay |18j , consists in approximating X using the Euler scheme with a uniformly 
spaced time grid. It suffers from two difficulties: first, for a general Levy measure 
V, there is no available algorithm to simulate the increments of the driving Levy 
process and second, a large jump of Z occurring between two discretization 
points can lead to a large discretization error. 

With the aim of resolving these problems, Rubenthaler [19] (see also Bruti- 
Liberati and Platen 0] and Mordecki et al [13] in the context of finite intensity 
Levy processes) introduced the idea of replacing the driving process Z by a suit- 
able compound Poisson approximation and placing the discretization points at 
the jump times of the compound Poisson process. This approach is problematic 
when the jump activity of the driving Levy process Z is strong, that is, the 
Levy measure has a strong singularity at zero. 

In Kohatsu-Tankov 0, the authors introduce and analyze a new approx- 
imation scheme in the case ct = 0, building on the ideas of Rubenthaler and 
Asmussen-Rosinski [2]. The idea is to replace the driving process Z by an ap- 
proximating process Z*^, which incorporates all jumps of Z bigger than e and 
approximates the jumps of Z smaller than e with a suitable chosen Brownian 
motion, matching the second moment of Z . The solution to the contiunuous 
SDE between the jump times can then be approximated with a suitable high 
order scheme. More recently, a similar approximation was used in the context 
of multilevel Monte Carlo schemes for Levy-driven SDKs [5] . 

Although the previous approach improves the rates of convergence obtained 
with Rubenthaler's scheme, there are limits on how well the small jumps of 
a Levy process can be approximated by a Brownian motion (think of non- 
symmetric Levy processes) . In Tankov |21) , the author presented a new scheme 
in the case u = based on approximating Z by a finite intensity Levy process, 
which incorporates all jumps bigger than e and matches a given number of 
moments of Z with an additional compound Poisson term. The main advantages 
of this approach are that the schemes are very easy to implement, because the 
driving process is piecewise deterministic, and that one can, in specific cases, 
obtain arbitrarily high order of convergence by matching a sufficiently large 
number of moments of Z. 

In this paper we are interested in two aspects of approximation schemes 
for Levy driven SDE's. First, in many of the previously mentioned schemes 
one assumes that there is no Brownian motion component in the equation ([I]) 
(i.e. cr = 0). The reason for this was that the speed of convergence of the ap- 
proximating scheme for the jump component is fast and therefore it was not 
clear how to match this speed with the approximation of the Brownian com- 
ponent without wasting computing resources. Furthermore the fact that the 
equation does not have a Brownian component facilitates the error analysis and 
the implementation of the scheme because the SDE between jumps is deter- 
ministic, as in [21], or can be treated as a deterministic equation perturbed by 
a small noise term as in [9]. On the other hand, recent developments in the 
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area of weak approximations for continuous SDE's [T51 [H] allow for high order 
approximations of the Brownian component. Therefore one may expect that 
the right combination of these approximation techniques with suitable jump 
adapted approximation schemes for pure jump SDE's can be achieved. 

Our second goal is a systematic study of the new moment-matching ap- 
proximation schemes introduced in |21| , with the objective of designing optimal 
compound Poisson approximations and studying their convergence in a more 
general setting. 

In this article, we show that the mathematical framework developed in 
Tanaka-Kohatsu [5D] is the appropriate tool in order to deal with the general 
situation (cr ^ 0). However, it needs to be adapted to the present setting where 
the partition is random while in [20], the partition is fixed. This framework is 
based on semigroup decompositions, which allow the study of a complex gener- 
ator by decomposing it into simple components. The error estimate is obtained 
by a local analysis of each component. 

In the resulting error bound, the contributions of the compound Poisson ap- 
proximation and of the discretization scheme for the Brownian part are separate 
and tractable. This allows to balance the two contributions by an appropriate 
choice of the order of the discretization scheme for the Brownian part, in order 
to minimize the computational time. On the other hand, this decomposition 
enables us to formulate the problem of choosing the compound Poisson ap- 
proximation as an optimization problem (minimizing the error bound). We 
characterize the optimal approximating process in the general case and provide 
explicit representation in specific situations. Often, the optimal solution is to 
keep all the jumps bigger than e and add an additional compound Poisson pro- 
cess to match the moment structure of the small jumps. Under a regularity 
assumption on the Levy measure, we show that this methodology can be used 
to construct approximations with arbitrarily high order of convergence. 

An interesting consequence of our analysis is that the Asmussen-Rosinski 
approach is not the optimal procedure to approximate the small jumps in the 
setting of weak convergence. We give a better procedure, which uses Levy 
measures with point masses to approximate the small jumps (see Remark 25 ) . 

In order to correctly describe the optimality aspect, let Xt be the unique 
solution of ([T]) but using Z as driving process instead of Z. Z is a finite activity 
Levy process with Levy measure D, which may have a Wiener component. Fur- 
thermore, let Xt be a computable approximation of Xt which shares the same 
jump times as X. The first objective is to find an upper bound for the difference 
Vi = E[f (Xi)] - E[f{Xi)] in terms oi X = D (R) < oo (the average number of 
partition intervals) and the moments of v — 9 and — v\. This part assumes 
then that the Brownian component can be simulated exactly. 

In the second part, we approximate the Brownian component and analyze 
the error Vi = E[/ (^i)] — E[/(Xi)]. To analyze 2?i, we extend the operator 
approach developed in |20j to jump-adapted random partitions. 

In conclusion, we find that we can express an upper bound for 2?i in terms of 
the moments oi u — D and \v—D\ and an upper bound for 2?i in terms of A. Now, 



3 



for fixed A (and, hence, Vi ) we consider P as a variable and minimize the upper 
bound for Pi , obtaining an optimal Levy measure P for the approximating finite 
intensity process Z. Once the optimal error is known as a function of A (this is 
done as a worse case analysis or in asymptotic form) one can identify the order 
of the approximation that is needed for the Brownian component. 

The paper is structured as follows. In Section [2] we introduce the notation. 
In Section [Sj we start introducing the assumptions in order to study the weak 
error of the approximations and we give the main error estimate, which will be 
the base for the study of optimal approximations. The expansion of the error 
is given in terms of A and the moments of v — D. 

The proof of the main error estimate is given in Sections |4.1| and |4.2[ which 
analyze, respectively, 2?i and Vi. In Section[5j we formulate the problem of find- 
ing the optimal compound Poisson approximation of Z as an optimization prob- 
lem, characterize its solution and prove an existence result. Explicit examples 
of solutions are given in Section |5.1[ and Section |5.2| analyzes the convergence 
rates of the resulting scheme. Specific algorithms and numerical illustrations 
are provided in Section [6j Finally, in the appendix we gather some technical 
lemmas. 

Throughout the article we use the Einstein notation of summation over dou- 
ble indices. Sy denotes the point mass measure at ?/ G K. Various positive 
constants are denoted by C or K with the dependence on various parameters. 
Their exact values may change from one line to the next without further men- 
tioning. 

2 Preliminaries and notation 

Let the process X = {Xt}t£[o^i] be the unique solution of the following d- 
dimensional SDE 

Xt=x+ f b (X,) ds+ f a (Xs) dBs + f h (Xs-) dZs, (2) 
Jo Jo Jo 

where b : R'^ ^ R'^ , h : R'^ ^ R'^ and a : R"^ ^ R'^'"' are {R'^) functions 
with bounded derivatives, B = {Bt}te[o,i] is a fc-dimcnsional standard Brownian 
motion and Z = {2^t}te[o,i] is a one dimensional Levy process (independent of 
B) with the following representation 

Zt^ f f yN {dy, ds)+ f f yN {dy, ds) , 

Jo J\y\<l JO J\y\>l 

N {dy, ds) = N (dy, ds) — v {dy) ds, 

where v is an infinite activity Levy measure, that is v {R) = +oo, and is a 
Poisson random measure on M x [0, oo) with intensity v {dy) x dt. 

Let X — {^t}tg[o,i] be the approximating process, which is the solution of 



4 



the SDE 

Xt=x^ f biXs)ds+ I a{X,)dB,+ f h{X,_)dZ,, (3) 
Jo Ja Jo 

where Z ~ {Zt}te[o,i] is a Levy process (independent of B) with the following 
representation 

Zt = flt + aWt+ f f yR{dy,ds)+ j f yN{dy,ds), 
Jo J\v\<i Jo J\v\>l 

N (dy, ds) = N {dy, ds) — v (dy) ds, 

where A = J^P (dy) < oo, ct^ > and iV is a Poisson random measure on 
M X [0,00) with intensity i> {dy) x ds and W = {M^f}te[o,i] is a standard k- 
dimensional Brownian motion independent of all the other processes. We assume 
that (/i, V, a) belongs to a set of possible approximation parameters denoted by 
A. Without loss of generality we may sometimes abuse the notation and write 
V € A to denote the Levy measure for which there exists p, and a so that 
V, a) S A. 

Note that, if we define 

b{x) = b{x) + h{x){p.- yv{dy)), 

J\v\<i 

then we can write 

Xt = x+ [ b{X,)ds+ [ a (X,) dB,+a [ h (X,) dW^ 
Jo Jo Jo 

h{X,_)yN{dy,ds). 



'0 JR 

Sometimes, the following flow notation will be useful 

Xt{s,x) = x + / b{Xu{s,x)) du+ / a {Xu{s,x)) dBu 

J S J S 

+ a [ h{Xu (s, x))dWu +11 h{Xu- {s, x))N {dy, ds) . 

J s J s Jk. 

Define the process 

Z{t,x) = x+[ b{%{t,x))du+ [ a{%{t,x))dB^ + a [ /i(?„(t,a;))rfW„ (4) 
Jt Jt Jt 

and the following operator 

{Ptf) {x)=E[m{0,x))]. 
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We consider the following stopping times 

= inf{t > f,_i : N (M, (T,_i,t]) ^0}, z G N, 

and the associated jump operators 

(57) (x) =E[/(x + ;i(a;) A%)], ieN 

Note that the stopping times Ti are well defined because A < oo and that 5' 
is independent of i because the jump sizes of a compound Poisson process are 
identically distributed. Still, we will keep this notation as it will help to keep 
track of the number of jumps. 

We will also assume that there exist a process X = {Xt]te[o,i] satisfying the 
following stochastic representation condition. 

Assumption 1 {STV) Assume that X satisfies 
for i £ N, where Pt is a linear operator. 

Remark 2 The process X and the linear operator Pt correspond to the scheme 
chosen to approximate the solution of equation ([3| between jumps. 

Recall that for each multi-index of order m, a = (ai, ad) G Z'^ we define 

\a\ := ai + ■ ■ ■ + ad — m. We also use the following notation /" — HiLi 

for any function / : M''' — > M.'^. We introduce the following spaces of functions. 

lid 



: the set of C" functions / 
a with < lal < m, 



such that for each multi-index 



fix) 



<C{a){l + \\x\n 



for some positive constant C (a) 



We will use the notation Cp :— Cp. In each C™ we consider the norm 



inf{C > 



5" 
dx"' 



fix) 



• Cft" : the set of C" functions / 
a with 1 < |a| < m, 

ga 

dx" 

for some positive constant C (a) 



< C(l-t-||a;f),0< |a| <m,xeR}. 
: M'' ^ M such that for each multi-index 



fix) 



<C{a) 
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Assumption 3 CHn) h,h,cF £ CJ^, J \yf" v (dy) < oo and sup^^^_^ J jyl^" (dy) < 
oo. 

Assumption 4 ('H'j) h,b,(T £ CJ^, J \y\'' v {dy) < oo andsup^^_^ J D (dy) < 
oo for all k > 1. 



In fact, all the results up to Section 4.2 only use moments up to power n + 1 



when we assume (T-Ln)- Still, in applications, in order for the continuous high- 
order scheme to satisfy the assumption {TZ{m,Sm)) (see below), the moments 
of order at least 2n are required. For this reason, we prefer this version of the 
assumptions. 

3 Weak error estimate 

Our next objective is to establish the main error estimate of this paper. In order 
to do this, we need to introduce a modification of the framework introduced in 



[20] in the next section. The error estimate will then be given in Section 3.2 



3.1 Framework for weak approximation of operator com- 
positions 

To simplify the notation, we define the non commutative product of operators 
as follows. Given a finite number of linear operators A^, A", we define 

n 

1=1 

Suppose we are given two sequences of linear operators {P^ }i>i and {Ql}i>i, 
t G [0, 1]. Furthermore, assume that for each i € N, Ql approximates in some 
sense to be defined later (see Assumption [7| . Given a partition tt = {0 = to < 
ti < • • • < tn-i < tn = 1}, we define its norm as |7r|„ := supj^j^ „{ti — 
Now, we would like to estimate the following quantity 

p,\Pl.t, ■ ■ ■ Pi-t,.^Ji^) - QlQl-t^ ■ ■ ■ Qi~u-^f{^) ■ 

In order to achieve this goal, we wih make use of the foUowing expansion 

n n 

n /k—1 n \ 

= Y.[X{Qu-U-APt-t.-.-Ql-t.-.) n Pl-U-Afi-)- (5) 

k=l \i=l i=k+l ) 

Hence, if we have a good norm estimates of 0^=1^ Pu-ti^i ^^'^ Yli=k+i Qu-ti^i 
then we can expect that n"=i Qu-ti^if i^) approximates well n"=i Pu-ti^if (^) • 

From now on, PI : Up>oCp -> Up>oCp, i £ N is a hnear operator for t £ [0, 1] 
and Ql : Up>oCp — > Up>oCp, i e N is a linear operator for t £ [0, 1]. 
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Assumption 5 (Ado) For all i gN, if f G Cp with p>2, then Q\f G Cp and 

sup ||0j/||^ <KiA)\\f\\c^, 
te[o,i] " 

for some constant K [A) > 0. Furthermore, we assume < Qlf {x) < Q\g {x) 
whenever < f < g and QJIr (x) = 1r (x). 

Assumption 6 {A4) For all i £ N, Ql satisfies {Ado) o,nd for each fp {x) := 

\xf (peN), 

Qlfp {x)<{l + K {A,p) t) fp (x) + K' {A,p) t 
for some positive constants K [A, p) and K' [A, p) . 

For m e N, : [0, 1] — )• R+ denotes an increasing function satisfying 

5m it) „ 

limsup q- = (J. 

Usually, we have Sm (t) = t™. 

Assumption 7 (7^(m,(5„)) For all i G N, define Errj = ErrJ"'' = - Q\. 
For each p > 2, there exists a constant q = q {m,p) such that if f G C™ with 
m* >2m + 2 then 

||Errj/|| <K{A,m)tSm{t)\\f\\cm' , 

for all t G [0, 1] . 

Assumption 8 {M p) If f G C™ one has that for k = 1, ...,n — 1 



sup 

{tk+i,-,t„)e[o,ir 



n pu 

!=fe+l 



<C{A) WfWcm. 

P 



Lemma 9 Under assumption {A4) , the operators {Ql}i>i satisfy 
sup sup max TT (31,_t,_j / (a;) < oo, 

A n i<k<n \ / 

for any positive function f G Cp, p > and |7r|„n < C for some positive 

constant C . 

Proof. Let fp (x) = \x\^^ for p e N. Using assumption [M.) , the monotonicity 
of the operators {Q\}i>i and that these operators are the identity on constants, 
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we have 
fc-i 



i=l 

<il + K iA,p) ih^i - tk-2)) (jl Ql-u_,^ fp i^) + K' {A,p) - tfc_2) 

<{l + K {A,p) kU (\iQl-u_}j fp {x) + K' {Ap) \An , 

with constants K {A,p) and K' {A^p) that do not depend on T:,x,k,n. Since 
(1 + K{A,p) |7r|„)''^^ < e'^-^(-^^P\ by induction follows that 

sup sup max TT Ql._t. \fix)<K' {A,p) e^(-^'f) (1 + \x\^P) < oo. 

_4 „ l<k<n \ -\^ ' ' I ^ ^ 



Theorem 10 Assume {A4) for PI and Q\ and {TZ (m, Sm)) ■ Then for any f € 
^2(m+i)^ i/iere exists a constant K — K {x^A^p) > such that 



k=l 



i=l i=l 

Proof. Let / € Cp'""^'^^\ Using the expansion ([5]), we have 

n n 

n /k—1 n 
k=l \i=l i=k+l 

Using assumption [TZ {m,Sm)) and {A4p) , we obtain 

/■ n \ 



< K |l/|lp2(™+i) ^ [tk - ifc-i) 5m {tk - tk-i) . 



< K {A, m) {tk - tk-i) Srn {tk - tk-i) (1 + |a;n 



i=k+l 



-^2(m + l) 



< K {A, m) {tk ~ tk-i) 5ra {tk - tk-i) (1 + |a;n ||/|L2(„+i) . 



Now, Lemma |9] yields 



i=k+l 



k-1 



< K (A m) (tk ~ h-i) 5,n {tk - tk-i) ||/|lc.^(™+i) n Qu-U-, ((1 + Nl')) 

i=l 

< K {x,A,m) {tk - tk-i) Sjn {tk - tk-i) ||/||p2(™+i) . 
Finally, adding up the estimates 



\{Pl-u.J{^)-\{Qu^u^J{^) 

n 

< K {x,A, m) |l/|1^2(™+i) ^ (tk - tk-i) 5jn {tk - tk-i) 



fe=i 



3.2 Main error estimate 

Theorem 11 Let X = {^t}tg[o,i] ^6 ^ process satisfying assumption {STZ) . 
Assume that the operators PI := S'^^^Pt and Q\ := S^^^Pt satisfy assumptions 
{M) and {TZ (m, Sm)) ,m>2. 

i) Assume {Un+i) and f G Cp^"+^' n C^"+\n > 2,p > 2. Then there exist 
positive constants K{x,A, m) and Ci{x), i = 1, n + 1 such that 

|E[/(Xi)]-E[/(Xi)]| 



< Ci {x) 



lal>i 



y{i' ~ v) (dy) ~ fi 



+ C2{x) 



y^{v ~ v) {dy) - a' 



+ / y^{v-v){dy) 

i=3 ^ 



+ Cn+i {x) I |yr+' - P| {dy) + K {x, A, m) ||/|L2(,„+i, A" 



ii) Assume {Hn+i) ^''^'^ f ^ 



(2(m+l))V(n+l) 



n > 2,p > 2. Then there exist 
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positive constants K(x,A, m) and Ci{x), i — 1, n + 1 such that 
|E[/(Xi)]-E[/(Xi)]| 



<C,ix)\\f\\^n 

P 

n 



y{v - v) {dy) - n 



|y|>i 



+ C2{x)\\f\\ 



y'^{v-v) {dy)-a^ 



i=3 



y\v~v) {dy) 



n+p+l 



{dy)} 



+ C„+i(x)|l/||^..i{/ \yr^\:.-D\{dy)+ / |y| 
+ K{x,A,m)\\f\\ A-™. 

Op 

Proof. Follows from Theorems [TSl and [TSl ■ 

Example 12 The first simple example of application of the above result is to 
parametrize the set A by a parameter e € (0, 1] so that: 



y{iy - Ve) {dy) 



\y\>i 



o^^al= [ y\u-y,){dy), 
Jr 

D{dy) EE v^{dy) = \i^\y\y^}v{dy). 



Take Pt = to be the operator associated with a one step Euler scheme, so 
that the overall approximation consists in applying the Euler scheme between 
the jumps of Z . Then the above result reads 

|E[/(Xi)]-E[/(Xf)]|<C3(a;) / \y\' u {dy) + K {x)\\f\\c. 

When (7 = 0, this result corresponds to Theorem 2 in fV^. 
In the particular case of an a-stabe-like Levy process with Levy density 
1^1 f+Q near zero, one obtains that the best convergence rate is X^^ for a < 1 

and the worse case is Ae for a — > 2. 

Note that we could have applied high order schemes for Wiener driven SDEs 
in order to improve the last term above to Aj™- 

Additional examples, algorithms, and numerical illustrations will be given 
in Section [6l 



4 Proof of the main error estimate 

4.1 Estimation of I?i = E[f (Xi)] - E[f (Xi)] 

Thoughout this section we will use the notation u {t, x) = E[f{Xi{t, x))]. Some 
auxiliary properties of this function u {t, x) are established in Lemma 37 
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Theorem 13 i) Assume (Hn+i) and f E C" ,71 > 2. Then, we have the 
following expansion 

E[f (Xi) - / (Xi)] = Bldt \ [ y{v- v) {dy) - fi I 
Jo [Av\>i ) 

(6) 
(7) 

i=3^0 JR JO 

where 



+ / Bfdtl / y'{,y-D){dy)^a^ 



J2 t Bldt [ y\v-v){dy)+ [' B^'dt, 

i=3 -^0 Js. Jo 



Bl := E 



:= E 



\a\=t 

E 

|Q|=ra+l 



1 a|a| _ _ ^ 



and 



xh" {Xt)y''+\u-D) {dy) 



\Bl\<C,{x), * = 
|Sr'|<C„+i(a;) / \y\^+'\v-v\{dy). 



(8) 



where the constants Ci {x) ,z = l,...,n+l, do not depend on v. 

ii) Assume anc? / e Cp^^,n > 2. T/ien we have that the expansion ^ 

also holds with \Bl\ < Ci (x) , i = 1, n, and 



B^^+^dt 



< 11/11 / \yr+'W-9\{dy) 

[Jr 

+ I \yr^+'\i^-i>\{dy)]. 



where the constants Ci {x) , i = 1, n + 1, do not depend on v. 

Proof. To simplify the notation we will give the proof in the case d = k — \. 
Note that E[/ {X{)\ = E[/ (Xi (0, x))\ = u (0, a;) and 



E[/ (^1) - / (^1)] = E[m (0, x)-u (1, ^1)]. 
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Applying Ito formula to u (l, Xi) and taking into account the equation satisfied 
by u {t,x) (see Lemma [37]) , we have 

E[u (0,x) -u(l,Xi)] 
"1 du 



= E 

+ E 

-E 



dx 



{t,Xt)h{Xt) \ / y{j^-D)idy)-fi\dt 



|y|>i 



du 



r^2 ,1^2^ 



, {t, Xt + h (Xt) y)^u {t, Xt) ~ ^ (t, Xt) h {Xt) y\{v-D) (dy) dt 



t,Xt)h' {Xt)dt 



Making a Taylor expansion of order n > 2, we obtain 



E 



du 



4=2 



{t, Xt + h (Xt) y)-u {t, Xt) - — {t, Xt) h (Xt) y\{'^-i^) {dy) dt 



dx' 



■u{t,Xt)h' (Xt) y\v-v) {dy)dt 



E 



^^7. {t,Xt + dyh {Xt)) —^d6 \ X h^+^ {Xt) y-^^\v - v) {dy) dt 



Hence, collecting terms, we have 



E 



du 
dx 



(tXt) h {Xt) 



dt ■ 



+ ^jj\v-v) {dy)-a^]E 



y{v - v) {dy) + / yi^ {dy) - /i , 

|y|>i 

21 dx^' 



'\y\<i 
ju{t,Xt) h' {Xt)dt 



pi 



1=3 

r-l 



Il-M^,xt)h^xt) 



E 



1 Qn+l 

dx 



dt ]j\^-^) {dy) 

(t, Xt + Xyh {Xt)) ^^^P^de^ X K'^+^ {Xt) f'+\v - v) {dy) 



dt 



and we obtain the expansion ([7| . Under the assumption ("Hn+i), using Lemmas 



34 and 37 



one obtains the first inequality in (|8|). Similarly, if we assume {T-LnJ\-i) 



using Lemmas 34 and 37 one obtains the second inequality in (|l 

4.2 Estimation of Vi = E[f (Xi)] - E[/(Xi)] 
Lemma 14 For z G N U {0}, one has that 

^i^{T,<l<T, + i}f{Xl)] = E[l^f.^-i^^f.^^yPi_fJ{Xf^ 
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Proof. Define W'^ := a{Xf^,Ti, ...,Ti+i),i eNU{0},j = Then, on the 

set {fi<l} 

r-1 



E 



b{X,{t,x))ds 



+ a{Xs{t,x))dBs + a h{Xs{t,x))dWs) 



t=Ti,x=Xj 



\t=Ti, 



where in the last equality we have used that Xg {t, x) satisfies the same SDE as 
Ys(t,x) on Ti < t < 1 < Ti+i. Now applying Lemma 36 and the definition of 
(Ptf) (x) we obtain the result. ■ 

Remark 15 Applying the previous lemma with i — and using that S'^ is the 
identity operator we obtain that 

E[l{i<^^j/(Xi)] = E[l{i<^^j5°A/ (x)]. 

Proposition 16 For i E N, the following equality holds. 

E[l{T.<i<f.+,}/ (^i)] = E[l{j,^<i<^,^^j5°%5i%_f^ • • • S^P,_fJ{x)]. 

Proof. Define G''^ a{Xf^_,fi, ...,fi+i),i e N,j = 1, By Lemma 
the definition of the operator 5' we have that 

E[l{T.<l<T.+,}/(^l)] 
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and 



{T,<l<T, + i}^l~Ti 



= E[1 

= E[l{^^<i<^^^^}E[Fi_^J(X^^_ + h{XT^^)AZr,)\g"]] 
= mm<i<f,+^}S'Pi-tf{x)\t=f,..=XfJ 
= E[l{2^^<i<j.^^jS'*Pi_j.^/(Xj.^_(Ti_i,Xf._ J)] 
= K[l{f,<i<f,+i}S^Pi~fJiYT,iTi-i,Xf._J)]. 
Where in the last equality we have used that 

/ h{X,{f,.i,Xf^ XyN{dy,ds)^0. 

Reasoning analogously to the proof of Lemma [llj one has that 



E 



1 { Ti < 1< Ti + 1 } "5" A - Ti / ( ^Ti ( - 1 1 -'^T, _ 1 ) ) 
E[l{T.<l<T. + ,}E[^^A-T./(%m-l,Xj,^_J)|H^ 
E [l{T.<l<T. + ,}E[5^A-tJ(l't, {h-i,x))]U^^f„u- 
E l{Ti<l<fi_|_i}-FTi-Ti_i'S"-Pl-Ti/(-^fi_i)] ■ 



= Xr, 
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Iterating this procedure the result follows. ■ 
Now we need the following technical result. 

Proposition 17 We have that 

oo i+1 



^ ^ E [l{f.<i<f.+,} (TkAl- r,_i)"+'] < C (m) A- 



i=0 k=l 

Proof. From Lemma 11 in [S', one has that 

r-l 



E 



{t ~ T] (t))"' dt 







< C(m) A" 



where rj (t) — sup{Ti : Ti < t} and C (to) is a constant that only depends on m. 
We can write 



E 



{t-r]{t)rdt 



^E l{T^<i<f.+,} / {t-v{t)rdt 

^=0 L -I 
oo i+1 /.TfcAl 

= E E E[im<i<f.+,} / it~v {t)r dti 

4=0 k=l JTk-i 

and the result follows by integration. ■ 

The main result of this section is the following. 

Theorem 18 Let {^t}te[o,i] &e the process defined in ^ and {^t}tg[o,i] o, 
process satisfying assumption {STZ) . If the operators PI :— S'^^^Pt and Q\ := 
S''~^Pt associated to these processes satisfy assumptions {A4) and {TZ{m,Sm)) 
with 6jn(t) — t™. Then for any f € Cp^'"^"'^' there exists a constant K = 
K {x,A,p) > such that 

E[/(Xi)] -E[/(Xi)]| <K{x,A,m) A-" 

Proof. We can write 



E[/(Xi)]-E[/(Xi)] -E 



E[im<i<T.+,}(/(^i)-/(^i)) 



By Proposition 16 and assumption {STZ), we have 



E 



E[im<i<f.+,}(/(^i)-/(^i)) 



.i=0 



EE 

i=0 



= EiE 



i=0 



i+1 

n 

fe=i 
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Then, by Theorem 10 we obtain that 
|e[/(Xi)]-E[/(Xi)]| 



OO 

i=0 



E 



i+1 



i+1 



Tk-i n ^n/\i-n-i 

k=l k=l 

oo i+1 

< K {x,A,m) \\f\\cl^^^+^)Y.Y.^[^{T.<i<T.+i} {TkM~ fk-i) 5m {fu M ~ Tk- 

i=a k=i 

Then the result follows by Proposition [T7| ■ 



/(^) 



5 Optimal approximation of Levy measures 



In this section, we discuss the optimization of the error bound in Theorem 11 
i) with respect to the choice of the approximating Levy process Z. We would 
like to choose the parameters p. and a and the Levy measure D in order to make 
the first four terms in the expansion small, that is, we concentrate on 



Ci{x) 



y(v - y) {dy) - ^ 



\y\>i 



+ C2[x) 



y^{v - v) {dy) - a' 



+ Y.'^i {x) f y\v - v) [dy) + Cn+i {x) f 

,- o JR JR 



i=3 

Our approach is to take 
M = 

'|y|>i 

so that the expansion becomes 



+ a,+,{x) I \yr+'W~i?\{dy) 



y{v — v) {dy) and cr = 



(9) 



J2c^{x) [ y\v~v){dy) + C„+i (x) / 

JR JR 



+ C„+i {x) / |yr+^ \v - v\ {dy) , 



(see Remark 25 for an alternative choice of a) 



Next, we choose the Levy measure v in the class of measures for which 
the first sum is equal to zero and then optimize over v in this class with fixed 
intensity A = S'(M) < cx) in order to make the last term as small as possible. We 
will denote by M. the set of all positive finite measures on M. The problem of 
finding the optimal approximating Levy measure then takes the following form. 

Problem 19 (fin, a) Letv be a Levy measure onM admitting the first n moments, 
where n > 2, and define = j^y^v{dy), 1 < k < n. For any D E M define 
the functional 



J{D) 



\yny-my)- 
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The problem Qn.A, n > 2, consists in finding 



fnfA) := min J M 



(10) 



under the constraints 

v{dy) = K and I y''D{dy) = m^, fc = 2, ...,n— 1, (11) 

: Jm 

where A > niinpgAf^ ^ where we set by convention minpgjv/i f'(]R) = 0. 

The computation of mmp^M„ for n > 2 is a classical problem, known 
as the Hamburger problem. A summary of known results on this problem is 
provided in Appendix A. 



Remark 20 In explicit examples of Section 5. 1 . and in the general treatment 
of Section 5.2. we shall see that for the solutions of Qn.A that we will find, the 
n+p+i _ I'lidy) appearing in Theorem 



term J^\y 
lower order as A 



oo than /jj - v\{dy). 



11 



a) will always be of a 
erefore, the convergence 



rates of our schemes will be the same under (Hn+i) ciiT'd under (jHn^i)- 
Proposition 21 The problem fin. a admits a solution. 



Proof. By Corollary 32 there exist at least one measure satisfying the con- 
straints ( [Tl] ). For n > 3, we define by the set of all such measures. For 
n = 2, we define by the set of all measures D € Ai satisfying D{R) = A and 
S^y'^v{dy) < C, where 



C 



x'^v{dx) 



It is clear that minimum in ( 10 1 is the same as the minimum over the set 
for any n > 2. 
Define 



Ka:={yeR: \y\ <a), a> 0. 
By Chebyshev's inequality we have that 

D{R\K,) ^ [ V {dy) < \ [ y^D {dy) , e A/^, 

J{\y\>a} O, 

which yields the tightness of M,l^. By Prokhorov's theorem, we have that the set 
is relatively sequentially compact but, as is closed (see e.g.. Chapter 
VII in Doob |6j), we also have that is sequentially compact. The set {J{P) : 
P £ -^n} bounded from below and, hence, it has an infimum, say (A). 
Then, by the basic properties of the the infimum, we can find a sequence of real 
numbers of the form {J {i'k)}k>i converging to fn(A). As is sequentially 
compact we can always find a sequence {i^fe, }z>i that converges weakly to some 
i?* g AI^. But {J (Pfe, )};>!, being a subsequence of the convergent sequence 
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{J {'^k)}k>i, must converge to £„(A). Hence, we only need to prove the lower 
semicontinuity of the functional J, that is, if Dk converges weakly to D then 
liminffc_>oo J [vk) > J {v) ■ 

Let D € M^. By the Hahn decomposition theorem, there exist disjoint 
measurable sets and S~ such that U S~ = M, u — v is nonnegative on 

and nonpositive on . The functional J(j^) can be alternatively written as 

J{D)= sup / \yrf{y){i,-v){dy), 

|y|"r(y)(i.-P)(dy), with r(y)-l5+ -Is-, 

where L°° is the space of bounded measurable functions endowed with the es- 
sential supremum norm. This implies that 

J{^) > sup / \y\^f{y){v-v){dy), (12) 

/eCo.||/||<i 

where Co is the space of continuous functions with compact support. 

Fix e > 0. By the monotone convergence theorem there exists A e (l,oo) 
such that ^ 

f \yrrmv-my)<e. 

J -A 

Since the measure /i :— |y|"(i^ — P) is a finite measure on M, both measures 
in its Jordan decomposition are also finite and hence inner regular (see e.g. 
V.16 in pj). Therefore, we can find two closed sets B'^ C 5*+ n {—A, A) and 
C 5^ n {—A, A) such that /i is positive on 5+, negative on B~ and At(M\ 
{B+ US")) < 2e. By Lusin's theorem, we can find an interpolation between 
\b+ and 1b-. That is, a function f G Cq with ||/|| < 1 such that f{x) = 1 for 
X e B+, f{x) = -1 for X e B^ and f{x) = for x ^ {-A, A) with 

fi{xeR; \f-lB-+lB+\{x)>e}<e. 

Therefore, finally 

■m - / |yr/(y)(^ - ^)(dy) < £ + |y|"(/*(y) - fivW - ^){dy) < 3e, 



which, together with (12) means that 



J{i?)^ sup / \yrfiy)ii^~D){dy), 

/6Co.||/||<l. 



because the choice of e was arbitrary. 
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For a sequence (vk) which converges weakly to z/, we have, for every f €L Cq 
with 11/11 < 1: 

\yrf{y){u - D){dy) = limmf ^ l2/r/(y)(^^ - i^k){dy) 

<hminf sup j \v\" f{y){v - Vk){dy) 

feC(u\\f\\<iJR 

= hminf J{vk)- 

k 

Now, taking the sup with respect to / in the left-hand side, we obtain the 
desired result. ■ 

The following result provides a characterization of the solutions of ^n,K^ 
which will be useful in finding explicit representations for small n. 



Proposition 22 The measure v is a solution o/ ( 10 1 if and only if it satisfies 



the constraints (111, and there exists a piecewise polynomial function P{y) = 
Oq + ^^^2 ^iV^ + such that P{y) > for all y e ]R, a function a : M [0, 1] 
and a positive measure t onM. such that 

v{dy) = v{dy)li^p^y)^2\v\^}+a{y)v{dy)li^p^y)=2\y\^} + {T{dy) + v{dy))li^P(y)=o}- 

(13) 

Remark 23 If the measure v is absolutely continuous with respect to Lebesgue's 



measure, the expression ( 13 1 simplifies to 



v(dy) = v{dy)\{p(y)^2\y\^} + r(d?;)l{p(j,)=o}- 

Moreover, in the case n = 2q, q € N, P (y) is a polynomial and the measure t 
may always be taken to be an atomic measure with at most q atoms (because a 
positive polynomial of degree n ^ 2q has at most q distinct roots). 



Proof. A measure D* which satisfies the constraints ( 11 ) is a solution of ( |10[ ) if 
and only if there exists a vector of Lagrange multipliers {po,P2, ■ ■ ■ ,Pn) such that 
u* minimizes the Lagrangian C{v,p) over all measures v ^ AA, and C{v* ^p) > 
—oo. The Lagrangian for this problem takes the form (dropping the terms which 
do not depend on D): 

« „ n—l 

C{D,p)^ / \y\-\y - y\{dy) + / ^^(dy) (po + ^ 

JVL ^^2 

Set P{y) ~ Po-^Y^^=2 PiV^ ~^\y\"' ■ 1^6* 2/0 G be such that P (yo) < 0, and consider 
the family of measures Va {dy) = aSy^, where a > 0. Then, for any (po, ■■•,Pn)j 



C-{va,p)= / \yY'vQ{dy) + \a-an\\yo\^ + a\po + y^piy:Q 
•^K\{yo} V i=2 
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where ao = i'{{yo}). For a > ao, we have that 



C{iya,p)= / \y\''Mdy)~ao\yo\'' + aP{yo) -oo. 

Therefore, necessarily P{y) > for all y S M. Now, as before, let the Jordan 
decomposition of D — i/ he given hy 9 — u = jj,^ — ii~ , where /i+ and /x" are 
supported on disjoint measurable sets. Then, 

C{i?,p)= [ P{y)f,+ {dy)+ [ {2\yr-P{y))f,-{dy) + C, 

JR JR 

where C denotes the terms which do not depend on /i+ and Then, it is 
clear that at optimum, 

• //+ should be equal to a measure with support {y : P{y) = 0}. Therefore 
in general, there will be no uniqueness. 

• fi- = Oon {y:2\y\^-P{y)>0}. 

• jjL~ = uq on {y : 2|t/|" — P{y) < 0}. This follows because /x+ and ^,~ are 
supported on disjoint measurable sets and ii~ < v. 

• /i~ satisfies v — [T' > on {y : 2|y|" — P(y) = 0}. 

Combining these observations, we complete the proof. ■ 

Example 24 Let n = 2q, q G N, and v he absolutely continuous. To find an 
optimal measure for the problem iln,A we can use the following procedure. Use 
the following parametrization for P (y) and r (dy) : 

P{y) = iy-aif ■■■{y-aqf , 

q 

i=l 

Solve the following system of n nonlinear equations for {ai}'l^i and {on}1^^ : 

Q Q 



/ (dy) + = A, 

•'{iy-air-{y-a^)^>2y^i} 

[ /i/(rfy) = Va,a^ k = 2, ...,2q - 1. 



'{(2/-ai)2--(y-a,)2>2y2<,} 

Obviously, in general, the solution to this system can only be approximated nu- 
merically and this does not seem an easy task. For n < A, the solutions are 
quite explicit; they are discussed in the following section. 
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To complete the analysis we need to quantify the dependence of the optimal 
value of the error f„ (A) on A when A tends to infinity. This is achieved in the 
following section for small values of n and in Section [5.2| for general n, under a 
regularity assumption on the Levy measure. 



5.1 Explicit examples for small values of n 

Throughout this section we assume that the measure v is absolutely continuous 
with respect to the Lebesgue measure. 



The case n = 2. We use the characterization of Proposition |22] (see also 
Remark 23 1. The function P{y) is necessarily of the form P{y) — aa + y'^ for 
some ao > (otherwise the infimum of the Lagrangian would be — oo), and 
therefore the optimal solution is given by 

i^e (dy) = l{y2y^}iy (dy) , 
where e — e{A) solves v{{y^ > s}) = A. The approximation error f2(A) is given 

by 

f2(A) = J{Pe{A)) - / y^Hdy), 
which can go to zero at an arbitrarily slow rate as A — >■ cxd. 



The case n — 3. The function P{y) is now of the form P{y) — ao + '^2y'^ + \y\'^ , 
and the positivity constraint implies that P{y) is necessarily of the form 

Piy)^{y + s)iy-2ef, y>0 
P{y) = -{y + 2ef{y~e), y < 0, 

or, in other words, P{y) — \y\^ — 3ey^ + Ae^, for some e > 0. It is now easy to 
see that an optimal solution is given by 

i^e (dy) = l{|y|>e}i^ (dy) + ai6-2s + 012626, 
where e — e(A) solves 



and 



ai + a2 = — ^ / y V [dy) . 



4£ 



{\y\<e} 



The approximation error satisfies f3(A) — o(A ^/^) as A — > 00, since 

f3(A) - / \y\^y{dx)+2e{K) f y^i^idx) < 3e(A) / y^iy^dx) = o(e(A)) 

J\y\<6(A) Ay\<eiA) J\y\<e{A) 
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and 



lim £(A)^A = lime^ / i/{dy)+lim\f y'^u{dx) <\\.m ( y'^v{dy)=Q. 

A-i-cx) £^0 J\y\>e ^-1-0 4 J\y\<e C^O J\y\<c 

However, the scheme with n = 4 achieves a better rate with the same compu- 
tational cost. 



The case n = 4. The function P{y) is now of the form P (y) = ao + a2y^ + 
asy^ + y'^ and from the positivity constraint we then deduce that 

P (y) = {y- ef{y + ef = y^ - 2y\^ + 

for some £ > 0. Analyzing the function 2y'^ — P (y) — y^ — £^ + 2e^2/^ it is easy to 
check that {2y^ — P (y) > 0} = {|y| > ^^/^72--l}. Hence, the optimal solution 
is of the form 

i?, {dy) = '^{dy)\^^^^^^/^^y + ai<5-£ + a^d,, 

where the constants ai and a2 are determined from the moment constraints 
and satisfy 



ai 



a2 



2£3 



{\v\<ey/V2-l} 

y'^v {dy) + e 



y^v {dy) +e y'^v {dy) ) , 

h\yV 



[/|<£\/V^-1} 



/ 



y^v {dy) 



2*^^ \J{\v\<eVV2-l} " ' ■ J{\y\<e^/V2-l}' 

and e = e(A) is found from the intensity constraint F {e) = A, where 



F{e) 



{\y\>£y/V2-l} 



v{dy) 



e J{\y\<es/V^r} 



y'^u {dy) . 



Note that F is strictly decreasing, continuous, and satisfies lim^^o F (e) = +oo 
and lim£-|-_|_oo F {s) = 0, which ensures the existence of a unique solution for 
F (e) = A. Also note that 



{\v\<ey/ V2-1} 



y^u {dy) 



< e 



< £ 



1 



/ 

J{\y 



|<e\/v^-l} 



{l!/l<eVV2-l} 

y'^v {dy) 



y'^v{dy) 



which ensures the non negativity of a\,a2- 

The worst case convergence rate can be estimated similarly to the case n = 3 
and satisfies Si{K) = o(A~^) as A — )■ oo. As we shall see in the next section, 
in the presence of a more detailed information about the explosion of the Levy 
measure at zero, this convergence rate can be refined. 
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Remark 25 



1. The calculations of this section make it clear that as far as weak approxi- 
mations are concerned, the Asmussen-Rosinski approach of approximating 
the small jumps of a Levy process with a Brownian motion is not nec- 
essarily the only answer. In fact, the case n = 3 studied above leads 
to an approximation which is asymptotically equivalent to the Asmussen- 
Rosinski method and the case n — A leads to a scheme which converges at 
a faster rate, for the same computational cost. 

2. Instead of taking a — 0, one may choose a which makes the second term 
in ^ equal to zero, which leads, for n > 3, to the following optimization 



problem for v: 



£:;(A) min J(z?) 

ueM 



under the constraints 

v[dy) = A and I y^v{dy) = nik, fc = 3, . . . , n — 1. 
Jr 

This problem assumes the use of the Asmussen-Rosinski approach to match 
the second moment of v. The analysis of this problem can be carried out 
using the same tools described above and leads to similar results. 



5.2 Convergence rates for regularly varying Levy mea- 
sures 

The notion of regular variation provides a convenient tool to study the conver- 
gence of our moment matching schemes even in the cases when n is large and 
an explicit solution of (il,i,A) is not available. We refer to [5] for background on 
regular variation. 

As usual, we denote by Ra the class of regularly varying functions with index 
a (at zero or at infinity depending on the context). The following assumption, 
which is satisfied by many parametric Levy models used in practice (stable, 
tempered stable/CGMY, normal inverse Gaussian, generalized hyperbolic etc.) 
may be used to quantify the rate of explosion of the Levy measure near zero. 

Assumption 26 There exists a e (0,2), positive constants c_|_ and c_ with 
C-f- + c^ > and a function g G R-a (o,t zero) such that the Levy measure v 
satisfies 

h'{{x,(x)) c^g{x) and ;/((— oo, — x)) ^ c_g(x) asx],0, (Ra)- 

Theorem 27 Let n be even and let the Levy measure v satisfy the assumption 
{Ra). Then there exists a function /(A) with f £ Ri^n/a as K oo such that 
the error bound £„ (A) defined by ( 10 1 satisfies 



cJ{K) < f„(A) < c/(A) 
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for all A sufficiently large, and for some constants c, c with < c <c < oo. The 
function f is given explicitly by /(A) = {g^{A))"A, where g^ is a generalized 
inverse of the function g appearing in Assumption {Ra)- 

Remark 28 

1. The regular variation implies that as A — >■ oo, the error goes to zero as 
A^^S times a slowly varying factor (such as logarithm). To compute the 
explicit convergence rate, the exact form of the regularly varying function 
g must be known. For example, if g{x) = then 

/(A) ^ CA'~^ 

for some strictly positive constant C. 

2. In the case n = A it can be shown using similar methods that £'n(A) 
C/(A) for some strictly positive constant C . 

Proof. Throughout the proof, we let g = ^ . To obtain an upper bound on the 
error, we construct a, possibly suboptimal, measure satisfying the constraints 
which attains the desired rate. Let e > 0, and define 

Ve{dy) = v{dy)l{\y\y^} + De{dy), (14) 

where Vi;[dy) is the solution (minimizer) of the moment problem 

A^ := min{zy(M) -.v^M, I y^v{dy) = to|. A; = 2, . . . 
where we define m| :— j^^y^^^-^y^v{dy). Then, 

fn(A,)< JK):= / y^\v-ve\{dy) 
Jr 

< f y",.{dy) + f y'^DMj) = 2 / 2/V(dy), (15) 

J{\y\<e} JR J{\y\<e} 

where 

A, -.^D.iR) = ( v{dy) + K,. 

By Proposition [30} 

Ae = inf{mQ : {"^j^^}^ > for some 

On the other hand, the matrix {?7T-f+j}f j=i is (nonnegative) positive definite, 
because it is a moment matrix of a measure. Therefore, by Sylvester's criterion 
we can write 

Ag = inf{r7io : det({TO^_|_j}f ^^q) > for some m\} 
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and also 



But 



Ae < inf{m§ : det({TOf+^-l,+,^i}«^.^o) > 0} 



det({r<+jli+j^i},f ,^.^0) = m^det({TO^+j}f^j-^i) +det({m^+^li+j>i}f^^-^o) 
and therefore 

A, < 



|det({m; 



''i+j-'-i+j>l}i,j=o)\ 



By integration by parts and Karamata's theorem (Theorem 1.5.11 in we 
show that 



lim 



^0 /(coo) Kt^^) P-"' 



for all p > a. 



(16) 



and so 



lim sup 



A, 



< 



det({^l,+,>i}?^^.^o) 



^io 4|>eKd^)" det({^}f,^.^i) 



The matrix {Aij)2j^i = ( is positive definite because 



dx. 



Therefore, det ^ > and there exits a constant C < 00 such that 



A,<C 



\z\>e 



v{dz) 



for £ sufficiently small. 

To sum up, we have found that there exist two positive constants Ci and C2 
such that for e sufficiently small, 



■J{\y\<e} 

Hdy) + A, < C2 



iy{dy) 



(17) 



A 



lal>e 



|y|>£ 

v{dy). 



Ial>e 



Let A(e) := j^^^^^vidy) and e(A) inf{e : A(e) < A}. This function satisfies 
A(£(A)) < A, and since A(£) e as e J, 0, by Theorem 1.5.12 in [3], we also 

get that £(A) e R—i/a as A — > 00. 

Now, for a given A, consider the measure (14) withe = e(A/C2), and possibly 
an additional atom at to satisfy the intensity constraint. This measure satisfies 
the constraints of Problem (r2„,A) and, by (17), has error bounded by 

Ci£"(A/C2);^ ^ CiC^'/"-'A£"(A), 
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so that the upper bound of the theorem holds with /(A) = Ae"(A) G Ri^n/a- 
To compute the lower bound, observe that 

^ri(A) > min J (i>) , 
and the explicit optimal solution for the problem in the right-hand side is given 

by 

where e and ^ € [0, 1] are such that J^y^^^ ^{dy) + £,i'[{\y\ = e}) = A (cf Propo- 
sition 22 1, which means that in particular e = £(A) introduced above. On the 
other hand, the error functional associated to this solution satisfies 

J{v,) = f \yrw - ^sl(dy) > / \yrHdy) - ^Ae"(A), 

Jr J\y\<e 

which proves the lower bound. ■ 



6 Description of the algorithm and numerical 
results 

According to Section [5j our approach to find an optimal approximation for the 
Levy measure starts by setting p, = yi^^ — ^)(dy) and ct = 0. Hence, the 

solution of equation ([s]) between jumps satisfies the following equation 



Yt{x)^x+ h{Y,{x))ds+ a{Ys{x))dBs, (18) 
Jo Jo 



b{x) = b{x) + "fh{x), 



where 



7 = / yv{dy) - / y^{dy). 

h\y\>^} h\y\>i} 

This implies that the drift term of the continuous part will depend on D through 
the parameter 7. Therefore, once we have fixed D the optimal approximation of 
the Levy measure v, we need to choose a weak approximation method to solve 



equation ( 18 1 . We will consider the following approaches: 



Weak Taylor approximations: These methods are based on the Ito- 



Taylor expansion of the solution of (18 1. This expansion is the stochastic 



analogue of the classical Taylor expansion, where the role of polynomials is 
played by multiple iterated stochastic integrals. Truncating the expansion 
at a certain degree of the iterated integrals we obtain an approximation 
method with global order of convergence related to that degree, see Propo- 
sition 5.11.1 in [8]. We will consider the weak Taylor approximations with 
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global order of convergence 1,2 and 3, which we will denote by WTl, WT2 
and WT3. Although the method is conceptually simple to understand, it 
presents some difficulties in the implementation as we need to sample from 
the joint law of multiple stochastic integrals of different orders. This makes 
the method less appealing from a practical point of view, especially when 
the driving Brownian motion is multi-dimensional. 

• Kusuoka-Lyons-Victoir methods: These methods are also based on 
stochastic Taylor expansions. The idea is to approximate the expecta- 
tion under the Wiener measure by the expectation under a probability 
measure supported on a finite number of paths of finite variation. By 
construction, the expectations of the iterated Stratonovich integrals, up 
to a certain degree, under this new measure match the expectations of 
the corresponding iterated integrals under the Wiener measure. Using 
the Stratonovich- Taylor formula one can deduce that the approximations 
obtained have global order of convergence depending on the degree of the 
iterated integrals taken into account, see [12]. In particular we will con- 
sider the approximation schemes of degree 3 and 5, denoted by KLV3 and 
KLV5, which give, respectively, global order of convergence 1 and 2. De- 
riving and implementing these methods is not straightforward, see [7j for 
an account on these issues. 

• Ninomiya-Victoir method: The Ninomiya-Victoir method can be seen 
as a stochastic splitting method. The idea is to find suitable small time 
approximations of the semigroup associated to the solution of equation 



( 18 ) . These approximations are written in terms of weighted products 
(compositions) of simpler semigroups associated to the so called coordi- 
nate processes and are deduced using formal Taylor expansions of the 
semigroups involved. The main difference with respect to the classical 
splitting methods is that, in the stochastic case, we need to find appropri- 
ate stochastic representations of the semigroups in order to implement the 
Monte Carlo method. These representations involve solving or approxi- 
mating ODEs with random coefficients. We will consider the algorithm 
given by Ninomiya and Victoir in [14) , which has global order of conver- 
gence 2. 

Having fixed an optimal Levy measure and a weak approximation scheme 
for the continuous part we can apply the following algorithm to obtain a sample 
of Xi. 

Algorithm to generate a weak approximation of Xi 

Requires: 

The initial condition x. 

The optimal Levy measure v. 

The weak approximation method (y) , to solve Yt iv) ,t & 
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Compute A = i^(M) and 7 = y^dy) - I{\y\>i} vHdv) 

Set Tiast = 0, Xnew = Xq 

Simulate the next jump time T ~ Exp(A) 

While (T < 1 - Tiast) do 

{ 

Compute Y^'^iXnew) 

Simulate A, a jump from the Poisson random measure 
with Levy measure D 

Set X^iew — (^ne-uj) h(Yrp (^Xnew))^ 

Set Tiast =T 

Simulate the next jump time T ^ Exp(A) 

} 

Compute Y^^^^^^{Xnetu) 
Set Xi — Y-^_rp^^^^^(yXnew) 

Return Xf"^ 



Applying, independently, the previous algorithm AI times we obtain a se- 
quence {xY^'^}i=i M and the Monte Carlo estimator of E[/(Xi)] is given 

by 

M 



M 

i=l 



We end this section with some numerical examples. We evaluate E[/(Xi)], 
where X is the solution of equation ([T]) with b{x) = ^oh{x) and a{x) = aoh{x). 
To approximate the Levy process, we use the optimal schemes presented in 
section [STT] with n = 2, n = 3 and n = 4, and denoted, respectively, by 0A2, 
0A3 and 0A4 in the examples below. For solving the continuous SDE between 
the times of jumps, we use the schemes WTl, WT2, WT3, KLV3, KLV5 and 
NV mentioned above. Finally, the process Z is taken to be a CGMY process, 
which is a Levy process with no diffusion component and Levy density of the 
form 

,e-^-l-ll,<o + e-^+l-ll,>o 



i^{x) ^ C- 



\l+a 



The third component of the characteristic triplet is chosen in such way that 
Z becomes a martingale. An algorithm for simulating the increments of Z is 
available [16 , which makes it possible to compare our methods to the traditional 
Euler scheme. Also, this process satisfies the assumption (Ra) of the previous 
section, and allows us to illustrate the dependence of the convergence rates on 
the parameter a. Actually, combining Theorems[TT]and[27]we have the following 
result. 



Theorem 29 Assume the hypotheses in Theorems 11 ii) and 21 , and choose 



7^ = and n — — i^) {dy) ■ Then, for n even, we have that there exist 
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positive constants K{x,A,m),C{x) and a slowly varying function I such that 

\E[f (Xi)]-E[/(Xi)]| 
< C(x)||/||p„+W(A)Ai-5 +if (a;,A™)ll/llc^(™+i)A-'", 

where A = 

We use 10^ simulation paths in all examples. For the Euler scheme, all values 
are computed using the same set of paths with 1, 2, 4, 8, 16, 32, 64, 128 and 256 
discretization intervals. For the optimal schemes, different paths are used for 
each point on the graph, and the different points are obtained by choosing the 
values of the parameter e which correspond to the values of := v{dx) 
in the range [0.5, 1, 2, 4, 8, 16, 32]. Also, the computing time for each point has 
been normalized by the standard deviation of the MC estimate, so that the 
times for all points correspond to the time required to get a standard deviation 
of 0.001. The variance of the MC estimate is about the same for all values 
computed with the optimal schemes. For the Euler scheme, the variance may 
be different, because, on one hand, the simulation method from |16j makes use 
of a probability change which increases variance, and on the other hand, we use 
a variance reduction techique for the Euler scheme (by taking E[/(a; + h{x)Zi)] 
as control variate) but not for the other schemes. In all the numerical examples 
below we take 70 = 0.5, ctq = 0.3, A+ — 3.5 and A_ — 2. Furthermore, for data 
set I, we take C = 0.5 and a = 0.5 (finite variation jumps) and for data set II 
we take C = 0.1 and a = 1.5 (infinite variation jumps). These two choices yield 
approximately the same variance of Xi and allow us to understand the effect of 
a on the convergence rate. 

For our first example, we take h{x) — x and j(x) = x. In this case, X 
is simply the stochastic exponential of 7oi + cr^Wt + -^t, and the exact value 
of E[/(Xi)] can be computed explicitly: E\l{Xx)\ = ■ Figure [l] plots the 
errors of the KLV schemes of different degrees and the NV scheme on a log-log 
scale for data sets I and II. In this case, the three approximations of the Levy 
measure, 0A2, 0A3 and 0A4, have very similar performance and we only plot 
the results for 0A2. This happens because with the choice j(x) = h{x) = x, we 
have E[/(Xi)] = E[/(Xi)] as soon as the approximation scheme for the Levy 
measure preserves the expectation of the Levy process, which is the case for 
all three approximation schemes OAl, 0A2 and OAS. In other words, for this 
choice of / and ft., the approximation of the Levy measure does not introduce 
any error. The error is therefore exclusively determined by the approximation 
scheme which is used between the jump times. However, in this case, the KLV 
and NV methods perfom so well that all the errors are below the statistical 
error due to the Monte Carlo method and it is not even possible to identify the 
actual order of convergence. 
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Figure 1: Errors of the cubature-based schemes for h(x) = x and f{x) = x. 
Left: parameters from data set I. Right: parameters from data set II. 

In our second example, we take h{x) = x still and f{x) — x^ . The exact 
value of E[/(Xi)] can also be computed explicitly and is now equal to 

E[X|] = E[£(2Z + [Z, Z])t] = exp{£;[2ZT] + E[[Z, Z]t]} 

= cxp |27or + <7^T + T y^Hdy)^ 

= exp {270T + a^T + TCV{2 - a)(A^-2 _^ y^-^)^ . 

Figure [2] plots the errors of the weak Taylor schemes of different orders on a 
log-log scale for data sets I and II, together with the theoretical error rates. 
In this case, one can clearly see the difference between the three schemes for 
approximating the Levy measure (0A2, OAS and 0A4) as well as the effect of 
the parameter a. 

For a = 0.5 (upper three graphs), the error of approximating the Levy 
measure is of order of A^~o = A~^ for 0A2, A~^ for OAS and A""^ for 0A4. 
Therefore, in these graphs, the global error is dominated by the one of approxi- 
mating the diffusion part: we observe a clear improvement going from WTl to 
WT2 and WTS, and no visible change going from 0A2 to OAS and 0A4. 

On the other hand, in the lower left graph, which corresponds to a = 1.5 and 
n — 2, the error of approximating the Levy measure is of order of A^^= — A~3 ^ 
which dominates the error of approximating the continuous SDE for any of the 
three weak Taylor schemes, and determines the slope of the curves in this graph. 
In this context, using the optimal scheme with n — S (lower middle graph) or 
n = 4 (lower right graph) leads to an substantial improvement of performance. 
In this case, we observe similar behavior for n = 3 and n — A because the 
Levy measure of Z is locally symmetric near zero, which means that 3-moment 
scheme and 4-moment scheme actually have the same convergence rate. 

The theoretical error rate of the Euler scheme is always ^ , which corresponds 
to the straight solid line on the graphs. The observed convergence rates appears 
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Figure 2: Errors of the weak Taylor schemes for h{x) = x and f{x) = x^ . Top: 
parameters from data set I. Bottom: parameters from data set II. 



slower than the theoretical prediction due to our variance reduction method, 
which has better performance when the number of discretization dates is small. 

A A moment matching problem 

In this section we present an auxiliary problem related with the moment match- 
ing of finite measures. 
We define 

Mn ■■={v eM: I y^v{dy) = m^, fc = 2, . . . , n}, 

where m^, k — 2, n are fixed real numbers. We want to compute infpg^^ D (M), 
i.e., the smallest intensity for which the moment constraints are feasible. This 
problem is very similar to the classical 'truncated Hamburger moment problem' 
and goes back to the works of Chebyshev, Markov and Stieltjes. The known 
results on an infinite interval can be summarized as follows |llj : 
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Proposition 30 Let n = 2q,q E N and let {mk}'^^Q be given . There exists 
a measure v £ A/„ with f'(M) = toq i] and only if the matrix 
nonnegative definite. 

Corollary 31 Let n — 2q, g G N, and let {wfclJJ^Q be given such that rrik = 
f^y^i'{dy), 2 < k < n for some nonnegative measure v. Then there exists a 
measure v e M„ with i?(]R) = mo if and only if det{{mi^j}'^ j^q) > 0. 



Proof. Using Proposition 30 it is enough to clieck tliat tlie tiie matrix {mi^j}1 



is nonnegative definite. ByTne definition of for k = 2, n we liave that the 
matrix {rui^j^l is nonnegative definite. Hence, by the Sylvester's criterion 
applied to the lower right corner minors of the matrix {Tii+j j^^-^Q, we have that 
in order for it to be nonnegative definite it is sufficient that det({mi+j}| ■ q) > 0. 



Corollary 32 For {rrikj^^^ ^■s Corollary 31, the set of values mo for which 
there exists a measure v € Mn with S'(K) = mo is of the form [mo,oo). 

The case when n is odd can be deduced from the previous one. 

Corollary 33 Let n = 2g + l,g e N. There exists a measure v e il/„ with 
p{M) = mo if and only if the matrix {mi+j}f^^Q is nonnegative definite for 
some mi £ M and m^+i G 

A simple matrix algebra computation then yields the following solutions for 
small n: 



n 


2 


3 


4 


5 


minpgTn,^ P (M) 








m4 
m4 


1714 



B Some useful lemmas on the solutions of SDEs 

In this section we will assume the notation established in the first section. 
Lemma 34 Assume that, for some p > 2, 

\yf V (dy) < oo, sup / \yf v (dy) < oo, 
i^eA Jm 

h,b,a G {M.).Then, there exists a constant C > 0, which does not depend on 
V, such that 



E 
E 



sup \Xtf 
.o<t<i 

sup At 
0<t<l 



<c(i + ixn, 
<c(i + ixn. 
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The proof of the this lemma is a standard generahzation of the proof for 
continuous sde's if one uses Kunita's second inequality (see Corollary 4.4.24 in 
Applebaum [1]). 

Lemma 35 Let p > 2 and for an integer n > 1 assume 

\y\"''v{dy) < oo, 

h,b,a ^ CJ^ (M). Then for any multi-index a with < \a\ < n we have 



E 



sup 

te[o,i] 



< oo. 



Proof. Follows from Theorem 70, Ch. V in p[3- ■ 

Using the time invariance of Levy processes one obtains the following result. 

Lemma 36 1. For <t < s < l,Xs (t,x) and Xs-t {0,x) have the same law. 

2. For < t < s < l,Ys(t,x) and Ys-t{0,x) have the same law, where 
Ys{t,x) is the process defined in Q. 

Lemma 37 Let u {t, x) = E[/ (Xi {t, x))] . 

(i) Assume (H„) and f e and bounded, with ?i > 2, . Then u e C^'" ([0, 1] x M) , 
are uniformly bounded for \a\ < n and u is a solution of the equation 



du 



+ {u{t,x + h (x) y)-u {t, x) - — (t, x) ht (x) y)v [dy) 

J\y\<l 



{u {t,x + h (x) y) ~ u (t, x)}v {dy) = 



(19a) 



'|y|>i 
u{l,x) = / (a;) 



(ii) Assume [U'^) and f e C?, with n > 2. Then u G C^^" ([0, 1] x M) , u is a 
solution of equation (19a) and there exists C < oo and p > with 



dx'^ 



(t,x) 



<c 11/11^. (i + i^n 



for all t e [0, 1], x e M and \a\ < 



Proof. The derivative satisfies 

ox 



du 
dxi 



(t,x) = E 



IL^ix,it,x))lxiit,x) 
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The interchange of the derivative and the expectations is justified using Lemma 
|35| Furthermore, one obtains by a direct estimation th e b oun ded ness under 
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and 



35 



The other 



{Tin) or the polynomial growth under using lemmas 
derivatives with respect to x are obtained by successive differentiation under the 
expectation and the derivative with respect to t is obtained from Ito's formula 



applied to /(Xi (t,x)) using Lemma 36 
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